bosaFS analyses
Pascale Caissy
2021-02-08
Section 1 Conditions environnementales régionales
Est-ce que les régions varient dans leurs conditions environnementales?
#Beaufort Sea: Mackenzie Shelf, Amundsen Gulf Maud
#Kitikmeot: Coronation Maud, Larsen Sound, Peel Sound
#Baffin Bay: North Water, Lancaster Sound, NW Baffin Bay
dfish$province <- with(dfish, ifelse(region %in% c("Mackenzie Shelf", "Amundsen Gulf Mouth"), "Beaufort Sea",
ifelse(region %in% c("Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound"), "Kitikmeot",
ifelse(region == "NEG", "Greenland Sea",
"NW Baffin Bay"))))
dfish$region_year <- as.factor(paste(dfish$region, dfish$year))
dfish$region_year <- ordered(dfish$region_year, levels = c("Mackenzie Shelf 2009","Mackenzie Shelf 2010","Mackenzie Shelf 2014","Mackenzie Shelf 2015",
"Amundsen Gulf Mouth 2015",
"Coronation Maud 2011" ,"Coronation Maud 2016","Coronation Maud 2017","Coronation Maud 2018",
"Larsen Sound - Victoria Strait 2010", "Larsen Sound - Victoria Strait 2016","Larsen Sound - Victoria Strait 2018",
"Peel Sound 2010","Peel Sound 2011","Peel Sound 2014","Peel Sound 2015","Peel Sound 2016",
"Lancaster Sound 2010","Lancaster Sound 2011","Lancaster Sound 2015" ,"Lancaster Sound 2016","Lancaster Sound 2017",
"North Water 2014","North Water 2016", "North Water 2018",
"West Baffin Bay 2015","West Baffin Bay 2016",
"NEG 2017"
))
station_data <- aggregate(data = dfish, unique_fish_id ~ region_year + station + surf_temp_degC + surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, length)
env_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, region_year, prof_mel) %>% PCA(quali.sup = c(4), graph = FALSE)
station_data$pc1 <- env_pca$ind$coord[, 1]
station_data$pc2 <- env_pca$ind$coord[, 2]
pca.vars <- env_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)
posArrowsX <- c(4, 4, 4,4)
posArrowsY<- c(4,4,4,4)
colours <- c(pals::stepped(n=24)[13:16], pals::polychrome(n=4)[4], pals::stepped(n=24)[5:11], pals::stepped2(n=20)[13], pals::stepped(n=24)[c(1:4,17:20)],pals::stepped3(n=16)[16],pals::stepped(n=24)[21:23], pals::stepped3(n=20)[c(5:6)], "red")
names(colours) <- levels(station_data$region_year)
gg_PCA_env <- ggplot() +
geom_point(data = station_data, aes(x = pc1, y = pc2, colour = region_year), size = 2.5)+
#scale_colour_manual(values = colours)+
coord_equal()+
theme_classic()+
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
legend.background = element_rect(fill=alpha('white', 0.4)))+
labs(x= paste("PC1 (", round(env_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
y= paste("PC2 (", round(env_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty = 2) +
geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
arrow = arrow(length = unit(0.025, "npc"), type = "open"),
lwd = 0.5) +
geom_text(data = pca.vars,
aes(x = Dim.1*posArrowsX, y = Dim.2*posArrowsY, label = vars),
check_overlap = FALSE)+
scale_color_manual(values = colours)
gg_PCA_env 
envzoo_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, NASC_zoo, region_year, prof_mel) %>% PCA(quali.sup = c(5), graph = FALSE)
station_data$pc1_zoo <- envzoo_pca$ind$coord[, 1]
station_data$pc2_zoo <- envzoo_pca$ind$coord[, 2]
pca.vars <- envzoo_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)
posArrowsX <- c(4, 4, 4,4,4)
posArrowsY<- c(4,4,4,4,4)
gg_PCA_envzoo <- ggplot() +
geom_point(data = station_data, aes(x = pc1_zoo, y = pc2_zoo, colour = region_year), size = 2.5)+
#scale_colour_manual(values = colours)+
coord_equal()+
theme_classic()+
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
legend.background = element_rect(fill=alpha('white', 0.4)))+
labs(x= paste("PC1 (", round(envzoo_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
y= paste("PC2 (", round(envzoo_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty = 2) +
geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
arrow = arrow(length = unit(0.025, "npc"), type = "open"),
lwd = 0.5) +
geom_text(data = pca.vars,
aes(x = Dim.1*posArrowsX, y = Dim.2*posArrowsY, label = vars),
check_overlap = FALSE)+
scale_color_manual(values = colours)
gg_PCA_envzoo
envfish_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, indstandard, region_year, prof_mel) %>% PCA(quali.sup = c(5), graph = FALSE)
station_data$pc1_fish <- envfish_pca$ind$coord[, 1]
station_data$pc2_fish <- envfish_pca$ind$coord[, 2]
pca.vars <- envfish_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)
posArrowsX <- c(4, 4, 4,4,4)
posArrowsY<- c(4,4,4,4,4)
gg_PCA_envfish <- ggplot() +
geom_point(data = station_data, aes(x = pc1_fish, y = pc2_fish, colour = region_year), size = 2.5)+
#scale_colour_manual(values = colours)+
coord_equal()+
theme_classic()+
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
legend.background = element_rect(fill=alpha('white', 0.4)))+
labs(x= paste("PC1 (", round(envfish_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
y= paste("PC2 (", round(envfish_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty = 2) +
geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
arrow = arrow(length = unit(0.025, "npc"), type = "open"),
lwd = 0.5) +
geom_text(data = pca.vars,
aes(x = Dim.1*posArrowsX, y = Dim.2*posArrowsY, label = vars),
check_overlap = FALSE)+
scale_color_manual(values = colours)
gg_PCA_envfish
envbio_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, NASC_zoo,indstandard, region_year, prof_mel) %>% PCA(quali.sup = c(6), graph = FALSE)
station_data$pc1_bio <- envbio_pca$ind$coord[, 1]
station_data$pc2_bio <- envbio_pca$ind$coord[, 2]
pca.vars <- envbio_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)
posArrowsX <- c(4, 4, 4,4,4,4)
posArrowsY<- c(4,4,4,4,4,4)
gg_PCA_bio <- ggplot() +
geom_point(data = station_data, aes(x = pc1_bio, y = pc2_bio, colour = region_year), size = 2.5)+
#scale_colour_manual(values = colours)+
coord_equal()+
theme_classic()+
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
legend.background = element_rect(fill=alpha('white', 0.4)))+
labs(x= paste("PC1 (", round(envbio_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
y= paste("PC2 (", round(envbio_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty = 2) +
geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
arrow = arrow(length = unit(0.025, "npc"), type = "open"),
lwd = 0.5) +
geom_text(data = pca.vars,
aes(x = Dim.1*posArrowsX, y = Dim.2*posArrowsY, label = vars),
check_overlap = FALSE)+
scale_color_manual(values = colours)
gg_PCA_bio
hist(aggregate(open_water_day ~ year + region, data = dfish, FUN = mean)$open_water_day)
plot(open_water_day ~ region, data = dfish, las = 2)
plot(open_water_day ~ as.factor(province), data = dfish)
aggregate(open_water_day ~ year + region +sampling_date + bu_date, data = dfish, FUN = mean)## year region sampling_date bu_date open_water_day
## 1 2009 Mackenzie Shelf 2009-08-05 2009-06-22 44
## 2 2009 Mackenzie Shelf 2009-08-15 2009-06-22 54
## 3 2010 Mackenzie Shelf 2010-08-16 2010-06-07 70
## 4 2010 Mackenzie Shelf 2010-08-19 2010-06-07 73
## 5 2010 Mackenzie Shelf 2010-08-21 2010-06-07 75
## 6 2010 Mackenzie Shelf 2010-08-22 2010-06-07 76
## 7 2010 Mackenzie Shelf 2010-08-23 2010-06-07 77
## 8 2010 Mackenzie Shelf 2010-08-24 2010-06-07 78
## 9 2010 Lancaster Sound 2010-08-07 2010-07-05 33
## 10 2010 Lancaster Sound 2010-08-08 2010-07-05 34
## 11 2010 Larsen Sound - Victoria Strait 2010-08-10 2010-08-09 1
## 12 2010 Peel Sound 2010-08-09 2010-08-16 -7
## 13 2011 Lancaster Sound 2011-08-04 2011-05-23 73
## 14 2011 Coronation Maud 2011-08-10 2011-07-18 23
## 15 2011 Peel Sound 2011-08-08 2011-07-25 14
## 16 2014 North Water 2014-08-01 2014-05-19 74
## 17 2014 North Water 2014-08-02 2014-05-19 75
## 18 2014 North Water 2014-08-06 2014-05-19 79
## 19 2014 Mackenzie Shelf 2014-08-23 2014-06-23 61
## 20 2014 Mackenzie Shelf 2014-08-24 2014-06-23 62
## 21 2014 Peel Sound 2014-08-10 2014-09-01 -22
## 22 2014 Peel Sound 2014-08-11 2014-09-01 -21
## 23 2015 Amundsen Gulf Mouth 2015-08-24 2015-05-11 105
## 24 2015 Lancaster Sound 2015-08-10 2015-05-18 84
## 25 2015 Lancaster Sound 2015-08-11 2015-05-18 85
## 26 2015 Lancaster Sound 2015-08-13 2015-05-18 87
## 27 2015 Mackenzie Shelf 2015-08-27 2015-05-25 94
## 28 2015 Mackenzie Shelf 2015-08-30 2015-05-25 97
## 29 2015 West Baffin Bay 2015-08-06 2015-07-06 31
## 30 2015 West Baffin Bay 2015-08-07 2015-07-06 32
## 31 2015 Peel Sound 2015-08-15 2015-08-10 5
## 32 2016 North Water 2016-08-04 2016-05-16 80
## 33 2016 North Water 2016-08-06 2016-05-16 82
## 34 2016 North Water 2016-08-07 2016-05-16 83
## 35 2016 North Water 2016-08-08 2016-05-16 84
## 36 2016 North Water 2016-08-09 2016-05-16 85
## 37 2016 North Water 2016-08-10 2016-05-16 86
## 38 2016 Lancaster Sound 2016-08-17 2016-05-23 86
## 39 2016 Lancaster Sound 2016-08-18 2016-05-23 87
## 40 2016 West Baffin Bay 2016-08-01 2016-07-04 28
## 41 2016 Coronation Maud 2016-08-21 2016-07-11 41
## 42 2016 Coronation Maud 2016-08-22 2016-07-11 42
## 43 2016 Coronation Maud 2016-08-23 2016-07-11 43
## 44 2016 Larsen Sound - Victoria Strait 2016-08-20 2016-07-25 26
## 45 2016 Peel Sound 2016-08-19 2016-08-01 18
## 46 2017 Lancaster Sound 2017-08-03 2017-05-22 73
## 47 2017 Lancaster Sound 2017-08-04 2017-05-22 74
## 48 2017 Lancaster Sound 2017-08-05 2017-05-22 75
## 49 2017 NEG 2017-08-25 2017-06-19 67
## 50 2017 NEG 2017-08-29 2017-06-19 71
## 51 2017 NEG 2017-08-31 2017-06-19 73
## 52 2017 NEG 2017-09-02 2017-06-19 75
## 53 2017 NEG 2017-09-04 2017-06-19 77
## 54 2017 NEG 2017-09-06 2017-06-19 79
## 55 2017 Coronation Maud 2017-08-08 2017-07-03 36
## 56 2017 Coronation Maud 2017-08-09 2017-07-03 37
## 57 2018 North Water 2018-08-28 2018-05-21 99
## 58 2018 North Water 2018-08-29 2018-05-21 100
## 59 2018 Coronation Maud 2018-08-21 2018-07-16 36
## 60 2018 Coronation Maud 2018-08-22 2018-07-16 37
## 61 2018 Larsen Sound - Victoria Strait 2018-08-19 2018-08-13 6lattice::xyplot(open_water_day ~ as.factor(year)|region, data = dfish)
hist(aggregate(surf_temp_degC ~ year + station, data = dfish, FUN = mean)$surf_temp_degC)
plot(surf_temp_degC ~ region, data = dfish, las = 2) #Beaucoup de variation Mackenzie Self
aggregate(surf_temp_degC ~ year + region +sampling_date + station + NASC_zoo + indstandard + surf_sal_kgm3 + open_water_day + bu_date, data = dfish, FUN = mean) %>% filter(region == "Mackenzie Shelf") ## year region sampling_date station NASC_zoo indstandard
## 1 2009 Mackenzie Shelf 2009-08-05 260 12.175651 0.266204537
## 2 2009 Mackenzie Shelf 2009-08-15 345 12.175651 0.010588621
## 3 2010 Mackenzie Shelf 2010-08-16 14 34.703792 0.011038416
## 4 2010 Mackenzie Shelf 2010-08-16 10 34.703792 0.022063208
## 5 2010 Mackenzie Shelf 2010-08-16 12 34.703792 0.013204152
## 6 2010 Mackenzie Shelf 2010-08-16 15 34.703792 0.005333019
## 7 2010 Mackenzie Shelf 2010-08-19 5 34.703792 0.007695008
## 8 2010 Mackenzie Shelf 2010-08-21 1 34.703792 0.038703415
## 9 2010 Mackenzie Shelf 2010-08-21 2 34.703792 0.016440352
## 10 2010 Mackenzie Shelf 2010-08-22 6 34.703792 0.006108735
## 11 2010 Mackenzie Shelf 2010-08-22 8 34.703792 0.007410152
## 12 2010 Mackenzie Shelf 2010-08-23 16 34.703792 0.012712653
## 13 2010 Mackenzie Shelf 2010-08-24 13 34.703792 0.010047687
## 14 2014 Mackenzie Shelf 2014-08-23 434 4.490901 0.005448571
## 15 2014 Mackenzie Shelf 2014-08-24 421 4.490901 0.001690748
## 16 2015 Mackenzie Shelf 2015-08-27 435 3.203247 0.002783745
## 17 2015 Mackenzie Shelf 2015-08-30 421 3.203247 0.000604846
## surf_sal_kgm3 open_water_day bu_date surf_temp_degC
## 1 29.46308 44 2009-06-22 4.1338462
## 2 27.75806 54 2009-06-22 2.1182560
## 3 26.13440 70 2010-06-07 9.2080000
## 4 27.16917 70 2010-06-07 8.6690000
## 5 27.94775 70 2010-06-07 8.0842500
## 6 28.05033 70 2010-06-07 8.0592222
## 7 25.85400 73 2010-06-07 9.7882000
## 8 25.92320 75 2010-06-07 10.1243000
## 9 27.73389 75 2010-06-07 8.0745789
## 10 25.05264 76 2010-06-07 9.9577374
## 11 27.83875 76 2010-06-07 9.2170000
## 12 28.38987 77 2010-06-07 6.9733871
## 13 26.32015 78 2010-06-07 9.8555385
## 14 30.72773 61 2014-06-23 5.2609500
## 15 26.73137 62 2014-06-23 0.1061250
## 16 28.34078 94 2015-05-25 3.9685500
## 17 27.09729 97 2015-05-25 -0.3764118plot(surf_temp_degC ~ as.factor(province), data = dfish)
lattice::xyplot(surf_temp_degC ~ as.factor(year)|region, data = dfish)
hist(aggregate(NASC_zoo ~ year+region, data = dfish, FUN = mean)$NASC_zoo)
plot(NASC_zoo ~ region, data = dfish, las = 2) #Beaucoup variation MS en 2010
plot(NASC_zoo ~ as.factor(province), data = dfish) #presque pas de zoo en mer du GL, et MS 2010 biaise vraiment la mer de Beaufort
plot(indstandard ~ region, data = dfish, las = 2) #Beaucoup variation LS -> c'est vraiment en 2010 qu'on a un pic
plot(log10(indstandard) ~ region, data = dfish, las = 2)
aggregate(indstandard ~ year + region +sampling_date + station, data = dfish, FUN = mean) %>% filter(region %in% c("Lancaster Sound")) ## year region sampling_date station indstandard
## 1 2010 Lancaster Sound 2010-08-07 301 0.90395209
## 2 2017 Lancaster Sound 2017-08-03 301 0.03870584
## 3 2017 Lancaster Sound 2017-08-04 304 0.01889668
## 4 2010 Lancaster Sound 2010-08-08 305 0.17822000
## 5 2016 Lancaster Sound 2016-08-18 305 0.02550317
## 6 2017 Lancaster Sound 2017-08-05 305 0.23590885
## 7 2011 Lancaster Sound 2011-08-04 332 0.09332530
## 8 2016 Lancaster Sound 2016-08-17 Allen Bay 0.03430469
## 9 2015 Lancaster Sound 2015-08-10 CAA-1 0.02836053
## 10 2015 Lancaster Sound 2015-08-10 CAA-2 0.11179698
## 11 2015 Lancaster Sound 2015-08-11 CAA-3 0.07652033
## 12 2015 Lancaster Sound 2015-08-13 CAA-4 0.10226513
## 13 2015 Lancaster Sound 2015-08-13 CAA-5 0.22170073plot(surf_sal_kgm3 ~ region, data = dfish, las = 2) #Beaucoup variation CM & PS, oz on a d'ailleurs des salinités assez faibles
lattice::xyplot(surf_sal_kgm3 ~ as.factor(year)|region, data = dfish)
aggregate(surf_sal_kgm3 ~ year + region +sampling_date + station, data = dfish, FUN = mean) %>% filter(region %in% c("Coronation Maud", "Peel Sound")) #2011 PS, 2017 et 2018 CM (QMG1-2-3-M) -> semble dépendre où était située la station QMG -> plus loin vers l'ouest = plus salé## year region sampling_date station surf_sal_kgm3
## 1 2014 Peel Sound 2014-08-10 309 30.46635
## 2 2010 Peel Sound 2010-08-09 310 28.32820
## 3 2011 Peel Sound 2011-08-08 310 19.93550
## 4 2014 Peel Sound 2014-08-11 310 30.28137
## 5 2016 Peel Sound 2016-08-19 310E 24.54000
## 6 2011 Coronation Maud 2011-08-10 314 25.85850
## 7 2016 Coronation Maud 2016-08-23 314 25.29824
## 8 2016 Coronation Maud 2016-08-23 316 27.13650
## 9 2015 Peel Sound 2015-08-15 CAA-7 29.13917
## 10 2016 Coronation Maud 2016-08-22 QMG-1 24.17040
## 11 2017 Coronation Maud 2017-08-09 QMG-1 21.01343
## 12 2018 Coronation Maud 2018-08-21 QMG-1 21.39480
## 13 2016 Coronation Maud 2016-08-21 QMG-2 20.82000
## 14 2017 Coronation Maud 2017-08-09 QMG-2 18.47175
## 15 2018 Coronation Maud 2018-08-21 QMG-2 23.47429
## 16 2016 Coronation Maud 2016-08-22 QMG-3 22.72938
## 17 2017 Coronation Maud 2017-08-09 QMG-3 27.54350
## 18 2016 Coronation Maud 2016-08-23 QMG-4 24.97867
## 19 2017 Coronation Maud 2017-08-09 QMG-4 27.97933
## 20 2018 Coronation Maud 2018-08-22 QMG-4 27.73485
## 21 2016 Coronation Maud 2016-08-21 QMG-M 19.80858
## 22 2017 Coronation Maud 2017-08-08 QMG-M 24.22665
## 23 2018 Coronation Maud 2018-08-22 QMG-M 25.57633plot(prof_mel ~ region, data = dfish, las = 2) #Beaucoup variation CM & PS, oz on a d'ailleurs des salinités assez faibles
range(dfish$prof_mel)## [1] 2.5 145.0lattice::xyplot(prof_mel ~ as.factor(year)|region, data = dfish)
Section 2: Facteurs reliés aux poissons eux-mêmes
Quel est l’éventail de taille de mes poissons et quel impact la taille semble-t’elle avoir?
La taille des poissons varie entre les régions, avec des larves particulièrement grosses dans l’Amundsen Gulf Mouth et Mackenzie Shelf ainsi qu’un énorme individu (61mm) dans Peel Sound. Les larves les plus petites sont en général celle de West Baffin Bay, Victoria Strait (avec un accroissement) et Peel Sound.
Les poissons échantillonnés dans une zone de débâcle hâtive sont généralement un peu plus gros, bien qu’il y ait encore des petits dans le pool. Ça reste que la majorité des poissons sont très petits lorsqu’ils sont collectés à une période concordant à ou près de la débâcle.
La quantité de carbone ingéré semble augmenter avec la taille. Par contre, comme on a beaucoup moins de grosses larves dans le pool, le signal est quand même dur à détecter (on a un gros biais pour les petites larves).
Quand on regarde le succès alimentaire en fonction de la taille, le nuage de point est assez confus. Généralement, le succès alimentaire brut est faible, peu importe la taille. Quand on transforme avec log10, on voit qu’il y a une légère augmentation en fonction de la taille. La relation linéaire entre ces deux variables est significative. Par contre, la variabilité semble à son plus haut lorsque le poisson a une très petite taille et une grande taille. Donc, ça indiquerait que parmi les très petits, on a dû en attraper qui seraient morts de famine bientôt OU qui avaient encore des réserves du sac vitellin et les gros seraient peut-être plus résistants à une famine variable, d’où la variabilité. Entre 12-13 et 25mm, le succès alimentaire augmente avec moins de variabilité. Donc, bien qu’on corrige FS en divisant par la masse, la taille demeure un facteur important de capture de proie.
#taille des poissons
plot(est_standard_length ~ region, data = dfish, las = 2)
tapply(dfish$est_standard_length, INDEX = list(dfish$region, dfish$year), FUN = mean)## 2009 2010 2011 2014 2015
## Amundsen Gulf Mouth NA NA NA NA 32.85720
## Coronation Maud NA NA 15.59700 NA NA
## Lancaster Sound NA 13.23085 13.62750 NA 16.59357
## Larsen Sound - Victoria Strait NA 11.30658 NA NA NA
## Mackenzie Shelf 12.6857 27.09488 NA 35.40000 31.54290
## NEG NA NA NA NA NA
## North Water NA NA NA 17.73251 NA
## Peel Sound NA 10.81823 17.42556 12.13674 11.36489
## West Baffin Bay NA NA NA NA 14.48783
## 2016 2017 2018
## Amundsen Gulf Mouth NA NA NA
## Coronation Maud 22.55682 16.26570 14.98442
## Lancaster Sound 17.96064 15.33450 NA
## Larsen Sound - Victoria Strait 14.74777 NA 11.03090
## Mackenzie Shelf NA NA NA
## NEG NA 19.12333 NA
## North Water 20.41569 NA 21.30636
## Peel Sound 12.85020 NA NA
## West Baffin Bay 11.19310 NA NAtapply(dfish$est_standard_length, INDEX = list(dfish$region, dfish$year), FUN = median)## 2009 2010 2011 2014 2015 2016 2017
## Amundsen Gulf Mouth NA NA NA NA 35.000 NA NA
## Coronation Maud NA NA 16.07 NA NA 18.906 16.1720
## Lancaster Sound NA 12.713 13.14 NA 15.857 18.438 14.6875
## Larsen Sound - Victoria Strait NA 11.431 NA NA NA 14.063 NA
## Mackenzie Shelf 12.307 27.259 NA 38.000 33.500 NA NA
## NEG NA NA NA NA NA NA 18.6500
## North Water NA NA NA 18.924 NA 21.997 NA
## Peel Sound NA 10.790 12.14 12.118 11.286 13.047 NA
## West Baffin Bay NA NA NA NA 13.857 11.136 NA
## 2018
## Amundsen Gulf Mouth NA
## Coronation Maud 15.1560
## Lancaster Sound NA
## Larsen Sound - Victoria Strait 10.7025
## Mackenzie Shelf NA
## NEG NA
## North Water 21.8750
## Peel Sound NA
## West Baffin Bay NAcor(dfish$est_standard_length, dfish$open_water_day)## [1] 0.4979243plot(dfish$est_standard_length~dfish$open_water_day)
lattice::xyplot(est_standard_length ~ open_water_day|region, data = dfish)
#taille et carbone
plot(dfish$carbon_mg ~ dfish$fish_WW)
plot(dfish$carbon_mg ~ dfish$est_standard_length)
#taille et succès alimentaire
plot((feeding_success)~est_standard_length, data = dfish)
plot(log10(feeding_success)~est_standard_length, data = dfish)
summary(lm(log10(feeding_success) ~ est_standard_length, data = dfish))##
## Call:
## lm(formula = log10(feeding_success) ~ est_standard_length, data = dfish)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.41983 -0.26364 0.06344 0.31037 1.15673
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.112901 0.058572 -53.147 < 2e-16 ***
## est_standard_length 0.009941 0.003028 3.283 0.00113 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4477 on 337 degrees of freedom
## Multiple R-squared: 0.031, Adjusted R-squared: 0.02812
## F-statistic: 10.78 on 1 and 337 DF, p-value: 0.001134plot(log10(feeding_success) ~ est_standard_length, data= dfish[dfish$est_standard_length>=25,]) #légère augmentation pour juvéniles en fonction de la condition, mais peu importe la condition chez les larves, le succès alimentaire reste assez variable (les résidus sont aussi moins dispersés)
plot(log10(feeding_success) ~ est_standard_length, data= dfish[dfish$est_standard_length<25,])
plot(log10(feeding_success) ~ fish_cond, data = dfish) #augmentation de FS
plot(log10(feeding_success) ~ fish_cond, data= dfish[dfish$est_standard_length>=25,])
plot(log10(feeding_success) ~ fish_cond, data= dfish[dfish$est_standard_length<25,]) 
plot(fish_cond~feeding_success, data = dfish)
plot(fish_cond~log10(feeding_success), data = dfish) # légère augmentation quand FS augmente
#pas vraiment de relation -> résidus pas assez dispersés (on n'a pas vraiment de valeurs extrêmes)
dfish$size_catego <- with(dfish, ifelse(est_standard_length<10, "<10",
ifelse(est_standard_length >=10 & est_standard_length<15, "10-15",
ifelse(est_standard_length>=15 & est_standard_length<20, "15-20",
ifelse(est_standard_length>=20 & est_standard_length<25, "20-25",
">25")))))
dfish$logFS <- log10(dfish$feeding_success)
xyplot(logFS ~ fish_cond|size_catego, data = dfish)
plot(fish_cond ~ est_standard_length, data = dfish) # la condition physique n'augmente pas vraiment avec la longueur du poisson
#FS et environnement
plot(logFS ~ NASC_zoo, data = dfish) #très légère augmentation mais rien de signifcatif
plot(logFS ~ surf_temp_degC, data = dfish) #légère augmentation
plot(logFS ~ surf_sal_kgm3, data = dfish) #pas grand chose
plot(logFS ~ open_water_day, data = dfish) #augmentation loglinéaire
dfish %>% select(logFS, open_water_day, surf_sal_kgm3, surf_temp_degC, NASC_zoo, fish_cond) %>% plot
#corrélation visible entre débâcle & température, débâcle et NASC_zoo et légère corrélation entre NASC_zoo et température
boxplot(logFS ~ region,data = dfish)
boxplot(logFS ~ province, data = dfish)
plot(feeding_success ~ region, data = dfish, las = 2) #pas mal similaire, mais toutes des valeurs bound à 0
plot(log10(feeding_success) ~ region, data = dfish, las = 2)
hist(dfish$feeding_success) #un beau neg binomial
hist(log10(dfish$feeding_success)) #normalish
lattice::bwplot(feeding_success ~ as.factor(year)|region, data = dfish)
lattice::bwplot(log10(feeding_success) ~ as.factor(year)|region, data = dfish) # variation LS, LV, MS, NW, PS
lattice::xyplot(log10(feeding_success) ~ (open_water_day)|region, data = dfish) #la variation individuelle est énorme toute région confondue
plot(feeding_success ~ as.factor(province), data = dfish)
plot(log10(feeding_success) ~ as.factor(province), data = dfish)
# Question 3 part 1
plot(prey_range_um/1000 ~ est_standard_length, data = dfish) # les grosses larves mangent à la fois de grosses et petites choses
plot(prey_min_um/1000 ~ est_standard_length, data = dfish)
plot(prey_max_um/1000 ~ est_standard_length, data = dfish)
plot(prey_median_um/1000 ~ est_standard_length, data = dfish) #vraiment rare que les larves, peu importe la taille, mangent de grosses affaires
plot(log10(feeding_success) ~ prey_median_um, data = dfish)
median(dfish[which(dfish$feeding_success > 0.001),"prey_median_um"])## [1] 292.5plot(log10(feeding_success) ~ prey_max_um, data = dfish)
plot(log10(feeding_success) ~ prey_min_um, data = dfish)
plot(log10(feeding_success) ~ prey_range_um, data = dfish)
plot((feeding_success) ~ prey_max_um, data = dfish)
plot((feeding_success) ~ prey_min_um, data = dfish)
plot((feeding_success) ~ prey_range_um, data = dfish)
#pas vraiment de relation entre taille des proies et FS
dfish$total_preys <- rowSums(dfish[,c(20:84)])
plot(total_preys ~ prey_max_um, data = dfish)
plot(total_preys ~ prey_min_um, data = dfish)
plot(total_preys~prey_range_um, data = dfish)
# Question 3 part 2
freqtable <- as.data.frame(table(gutcontent2$prey_id_comb,gutcontent2$region))
freqtable$Var2 <- ordered(freqtable$Var2, c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "NEG"))
color = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)]
colours_bar <- sample(color, size =15)
ggplot() +
geom_bar(data =freqtable, mapping = aes(fill=Var1, y=Freq, x=Var2), position="fill", stat="identity") +
scale_fill_manual(values =colours_bar)
freqtable2 <- as.data.frame(table(gutcontent2$prey_id_comb,gutcontent2$station))
ggplot() +
geom_bar(data =freqtable2, mapping = aes(fill=Var1, y=Freq, x=Var2), position="fill", stat="identity") +
scale_fill_manual(values =colours_bar)
carbon_input <- aggregate(carbon_mg ~ prey_id_comb + region, data = gutcontent2, FUN = sum)
ggplot() +
geom_bar(data =carbon_input, mapping = aes(fill=prey_id_comb, y=carbon_mg, x=region), position="fill", stat="identity") +
scale_fill_manual(values =colours_bar) +
labs(x = "", y = "")
carbon_input2 <- aggregate(carbon_mg ~ prey_id_comb + station, data = gutcontent2, FUN = sum)
ggplot() +
geom_bar(data =carbon_input2, mapping = aes(fill=prey_id_comb, y=carbon_mg, x=station), position="fill", stat="identity") +
scale_fill_manual(values =colours_bar) 
#en gros, bouffer des calanus = payant à souhaitPCA
fish_pca <- dpreys_comb %>% select(-unique_fish_id,-fish_WW, -est_height_anus,-carbon_mg, -surf_sal_kgm3, -station, -sampling_date, -year, -bu_date, -c( 17:18, 20:42)) %>%PCA(quali.sup = c(2), quanti.sup = c(7))
dpreys_comb$pc1 <- fish_pca$ind$coord[, 1]
dpreys_comb$pc2 <- fish_pca$ind$coord[, 2]
pca.vars <- fish_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)
pca.vars.m <- reshape2::melt(pca.vars, id.vars = "vars")
pca.var.sup <- fish_pca$quanti.sup$coord %>% data.frame
colours <- c("#6BCA54", "#FF5733", "#73C6B6", "#5E82D6", "#424949", "#A93226", "#C39BD3", "#FFC300", "#BDC3C7")
posArrowsX <- c(4, 4, 4, 3.25,4,4)
posArrowsY<- c(4.25, 3.5, 4.75, 4.25,4,4)
p_transf <- ggplot() +
geom_jitter(data = dpreys_comb, aes(x = pc1, y = pc2, colour = region), size = 2.5)+
coord_equal()+
theme_classic()+
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
legend.position = c(0.25,0.19),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size = 20),
legend.text = element_text(size = 22),
axis.title = element_text(size = 25),
legend.title = element_text(size =22),
legend.background = element_rect(fill=alpha('white', 0.4)))+
labs(x= paste("PC1 (", round(fish_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
y= paste("PC2 (", round(fish_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
geom_hline(yintercept = 0, lty = 2) +
geom_vline(xintercept = 0, lty = 2) +
geom_point(size = 2) +
geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*4, y = 0, yend = Dim.2*4),
arrow = arrow(length = unit(0.025, "npc"), type = "open"),
lwd = 0.5) +
geom_text(data = pca.vars,
aes(x = Dim.1*posArrowsX, y = Dim.2*posArrowsY,
label = c("LS","Jours depuis\ndébâcle","Abondance du\nzooplankton", "Température de\nsurface","Compé", "Condition physique")),
check_overlap = FALSE, size = 8.5)+
scale_color_manual(values = colours, name = "Region") +
geom_segment(data = pca.var.sup, aes(x = 0, xend = Dim.1*5, y = 0, yend = Dim.2*5),
arrow = arrow(length = unit(0.025, "npc"), type = "open"),
lwd = 0.5, lty = "dotted") +
geom_text(data = pca.var.sup,
aes(x = Dim.1*4, y = Dim.2*4,
label = "Succès\nalimentaire"), size = 8.5, col = '#868B8D')
p_transf
ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/PCA_fish.png", plot = p_transf, dpi = 500, width = 14, height = 10, bg = "transparent")CA
CA.preys <- dpreys_comb %>% select(22:42) %>% CA()
summary(CA.preys) #first two axes explains 24%, not really good
plot(CA.preys)
# PCA.prey <- dpreys_comb %>% select(-unique_fish_id, -est_height_anus, -fish_WW, -carbon_mg) %>% PCA(quanti.sup = c(7), quali.sup = c(2))
env <- dpreys_comb %>% select(-region, -unique_fish_id, -est_height_anus, -fish_WW, -carbon_mg, -surf_sal_kgm3, -est_standard_length, -c(13:34))
preys <- dpreys_comb %>% select(13:34) %>% decostand(method = "hellinger")
rda.preys <- rda(X = preys, Y = env, scale = TRUE)
plot(rda.preys)
summary(rda.preys) #not really good, two first axes explain only 11%...
preys_pres_abs <-dpreys_comb %>% select(13:34) %>% sapply(FUN = function(x)
ifelse(x > 0, 1, 0))
rda.preys.presabs <- rda(X = preys_pres_abs, Y = env, scale = TRUE)
plot(rda.preys.presabs)
summary(rda.preys.presabs)
PCA.prey <- dpreys_comb %>% select(-unique_fish_id, -est_height_anus, -fish_WW, -carbon_mg, -est_standard_length) %>% PCA(quanti.sup = c(6), quali.sup = c(1))Bray Curtis, MarineCusa, and others…….
wo data sets were generated for the multivariate analysis: Dataset A, where the gravimetric contribution of BDC in percent (%W) was calculated for each stomach by dividing prey taxa weight by total prey weight for each stomach and multiplying the result by 100;
and dataset B, in which the average BDC weight was calculated for each trawl by adding total prey taxa weight for each trawl, dividing total prey taxa weight by total prey weight for each trawl, and multiplying the result by 100. Dataset B calculation were performed separately for each polar cod size category
Sélection de modèles
#vérification ajustement des modèles lm----
lm.fs.untranf <- lm(feeding_success ~ open_water_day + fish_cond, data = dfish, na.action = na.omit)
plot(lm.fs.untranf)
lm.glob1 <- lm(log10(feeding_success) ~ open_water_day + fish_cond, data = dfish, na.action = na.omit)
summary(lm.glob1)##
## Call:
## lm(formula = log10(feeding_success) ~ open_water_day + fish_cond,
## data = dfish, na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.45771 -0.25781 0.05526 0.28660 1.24458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.115499 0.041676 -74.755 < 2e-16 ***
## open_water_day 0.003567 0.000692 5.155 4.33e-07 ***
## fish_cond 0.183561 0.073807 2.487 0.0134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4321 on 336 degrees of freedom
## Multiple R-squared: 0.1003, Adjusted R-squared: 0.09497
## F-statistic: 18.74 on 2 and 336 DF, p-value: 1.932e-08#pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/glob1_res.pdf")
plot(rstudent(lm.glob1) ~ fitted(lm.glob1), ylab = "Résidus de Student", xlab = "Valeurs prédites")
#dev.off()
plot(hatvalues(lm.glob1),
ylab = "Hat values",
xlab = "Observations", main = "Effet de levier")
abline(h = 2*3/324, lty = 2) #on a un effet de levier pour certaines valeurs, mais je ne crois pas qu'il soit trop
# pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/glob1_cook.pdf")
plot(cooks.distance(lm.glob1),
ylab = "Distance de Cook",
xlab = "Observations",
main = "Influence des observations")
abline(h = 3/(324-4), lty = 2) #rien en haut de 1, mais un bon nombre sont en haut du seuil
# dev.off()
lm.glob2 <- lm(log10(feeding_success) ~ surf_temp_degC + NASC_zoo + fish_cond, data = dfish, na.action = na.omit)
summary(lm.glob2)##
## Call:
## lm(formula = log10(feeding_success) ~ surf_temp_degC + NASC_zoo +
## fish_cond, data = dfish, na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4832 -0.2482 0.0573 0.3056 1.1629
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.9486490 0.0407859 -72.296 < 2e-16 ***
## surf_temp_degC 0.0024438 0.0132247 0.185 0.85350
## NASC_zoo 0.0004823 0.0039241 0.123 0.90226
## fish_cond 0.2391899 0.0791796 3.021 0.00271 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4494 on 335 degrees of freedom
## Multiple R-squared: 0.02951, Adjusted R-squared: 0.02082
## F-statistic: 3.395 on 3 and 335 DF, p-value: 0.01815# pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/glob2_res.pdf")
plot(rstudent(lm.glob2) ~ fitted(lm.glob2), ylab = "Résidus de Student", xlab = "Valeurs prédites")
# dev.off()
plot(hatvalues(lm.glob2), ylab = "Hat values",
xlab = "Observations", main = "Effet de levier")
abline(h = 2*4/324, lty = 2) #on a un effet de levier pour certaines valeurs, mais je ne crois pas qu'il soit trop
# pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/glob2_cook.pdf")
plot(cooks.distance(lm.glob2), ylab = "Distance de Cook", xlab = "Observations", main = "Influence des observations")
abline(h = 4/(324-4), lty = 2) #rien en haut de 1, mais un bon nombre sont en haut du seuil
dev.off()## null device
## 1#modèles mixtes----
lme.null <- lme(log10(feeding_success) ~1, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.null)## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 412.5423 424.0203 -203.2712
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1415781 0.4318548
##
## Fixed effects: log10(feeding_success) ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) -2.913366 0.05418298 330 -53.76903 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2166700 -0.5841590 0.1201451 0.6732820 3.0226108
##
## Number of Observations: 339
## Number of Groups: 9lme.glob1 <- lme(log10(feeding_success) ~ surf_temp_degC + NASC_zoo + fish_cond + indstandard, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.glob1) #serait mieux d'utiliser une transformation log10## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 413.8603 440.6423 -199.9302
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1347363 0.4279663
##
## Fixed effects: log10(feeding_success) ~ surf_temp_degC + NASC_zoo + fish_cond + indstandard
## Value Std.Error DF t-value p-value
## (Intercept) -2.8890499 0.06751992 326 -42.78811 0.0000
## surf_temp_degC 0.0008893 0.01568588 326 0.05670 0.9548
## NASC_zoo -0.0036067 0.00477970 326 -0.75458 0.4510
## fish_cond 0.1759375 0.07905116 326 2.22562 0.0267
## indstandard -0.1204569 0.22282177 326 -0.54060 0.5892
## Correlation:
## (Intr) srf__C NASC_z fsh_cn
## surf_temp_degC -0.480
## NASC_zoo -0.017 -0.588
## fish_cond 0.044 -0.179 0.191
## indstandard -0.276 0.132 0.007 -0.021
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.19707455 -0.51683825 0.09907741 0.64911501 3.07875340
##
## Number of Observations: 339
## Number of Groups: 9lme.glob2 <- lme(log10(feeding_success) ~ open_water_day + fish_cond, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.glob2)## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 392.4367 411.5667 -191.2183
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1226326 0.4176763
##
## Fixed effects: log10(feeding_success) ~ open_water_day + fish_cond
## Value Std.Error DF t-value p-value
## (Intercept) -3.165830 0.07516693 328 -42.11732 0.0000
## open_water_day 0.004910 0.00113264 328 4.33509 0.0000
## fish_cond 0.153737 0.07540482 328 2.03882 0.0423
## Correlation:
## (Intr) opn_w_
## open_water_day -0.766
## fish_cond 0.070 -0.097
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.1276007 -0.5938883 0.1068560 0.6323493 3.0566429
##
## Number of Observations: 339
## Number of Groups: 9performance::r2(lme.glob2)## # R2 for Mixed Models
##
## Conditional R2: 0.217
## Marginal R2: 0.150lme.glob3 <- lme(log10(feeding_success) ~ open_water_day + fish_cond +surf_temp_degC + NASC_zoo + indstandard, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.glob3)## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 392.1609 422.7689 -188.0805
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1486085 0.4121206
##
## Fixed effects: log10(feeding_success) ~ open_water_day + fish_cond + surf_temp_degC + NASC_zoo + indstandard
## Value Std.Error DF t-value p-value
## (Intercept) -3.182257 0.09219465 325 -34.51673 0.0000
## open_water_day 0.006638 0.00134656 325 4.92936 0.0000
## fish_cond 0.136769 0.07680118 325 1.78082 0.0759
## surf_temp_degC -0.013225 0.01558379 325 -0.84865 0.3967
## NASC_zoo -0.006531 0.00471398 325 -1.38551 0.1668
## indstandard 0.174968 0.22380914 325 0.78177 0.4349
## Correlation:
## (Intr) opn_w_ fsh_cn srf__C NASC_z
## open_water_day -0.653
## fish_cond 0.090 -0.096
## surf_temp_degC -0.217 -0.183 -0.149
## NASC_zoo 0.063 -0.114 0.194 -0.553
## indstandard -0.358 0.264 -0.045 0.067 -0.015
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.21435875 -0.60432496 0.08542066 0.64163558 3.16442861
##
## Number of Observations: 339
## Number of Groups: 9performance::r2(lme.glob3)## # R2 for Mixed Models
##
## Conditional R2: 0.286
## Marginal R2: 0.193lme.bu <- lme(log10(feeding_success) ~ open_water_day, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.bu)## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 394.5949 409.8989 -193.2975
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1280725 0.4199369
##
## Fixed effects: log10(feeding_success) ~ open_water_day
## Value Std.Error DF t-value p-value
## (Intercept) -3.179778 0.07694985 329 -41.32274 0
## open_water_day 0.005202 0.00114956 329 4.52555 0
## Correlation:
## (Intr)
## open_water_day -0.761
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.16166520 -0.57445494 0.09310238 0.66866541 2.98507700
##
## Number of Observations: 339
## Number of Groups: 9lme.cond <- lme(log10(feeding_success) ~ fish_cond, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.cond)## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 408.8741 424.1781 -200.4371
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1305161 0.4288933
##
## Fixed effects: log10(feeding_success) ~ fish_cond
## Value Std.Error DF t-value p-value
## (Intercept) -2.9155946 0.05090030 329 -57.28050 0.0000
## fish_cond 0.1845831 0.07704973 329 2.39564 0.0171
## Correlation:
## (Intr)
## fish_cond -0.006
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.16832072 -0.54686777 0.09007215 0.65614797 3.08828777
##
## Number of Observations: 339
## Number of Groups: 9lme.comp <- lme(log10(feeding_success) ~ indstandard, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.comp)## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 414.3573 429.6613 -203.1786
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1399442 0.4318376
##
## Fixed effects: log10(feeding_success) ~ indstandard
## Value Std.Error DF t-value p-value
## (Intercept) -2.9084324 0.05509209 329 -52.79220 0.0000
## indstandard -0.0951137 0.22091146 329 -0.43055 0.6671
## Correlation:
## (Intr)
## indstandard -0.217
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2231890 -0.5886669 0.1248675 0.6830543 2.9948472
##
## Number of Observations: 339
## Number of Groups: 9lme.NASC <- lme(log10(feeding_success) ~ NASC_zoo, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.NASC)## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 413.3037 428.6077 -202.6519
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1474387 0.4306864
##
## Fixed effects: log10(feeding_success) ~ NASC_zoo
## Value Std.Error DF t-value p-value
## (Intercept) -2.8862064 0.06073509 329 -47.52124 0.0000
## NASC_zoo -0.0043282 0.00386754 329 -1.11910 0.2639
## Correlation:
## (Intr)
## NASC_zoo -0.388
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.24785165 -0.57784643 0.08928769 0.65987904 3.05791490
##
## Number of Observations: 339
## Number of Groups: 9lme.temp <- lme(log10(feeding_success) ~ surf_temp_degC, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.temp)## Linear mixed-effects model fit by maximum likelihood
## Data: dfish
## AIC BIC logLik
## 414.4939 429.7979 -203.247
##
## Random effects:
## Formula: ~1 | region
## (Intercept) Residual
## StdDev: 0.1430711 0.4317298
##
## Fixed effects: log10(feeding_success) ~ surf_temp_degC
## Value Std.Error DF t-value p-value
## (Intercept) -2.9045914 0.06695295 329 -43.38258 0.0000
## surf_temp_degC -0.0027957 0.01259061 329 -0.22205 0.8244
## Correlation:
## (Intr)
## surf_temp_degC -0.577
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2347463 -0.5854568 0.1171700 0.6707042 3.0315922
##
## Number of Observations: 339
## Number of Groups: 9list.models <- list(null = lme.null, tempNASCFishcond = lme.glob1, buFishcond= lme.glob2, buOnly = lme.bu, tempONLY = lme.temp, condOnly = lme.cond, NASConly = lme.NASC, compONLY = lme.comp)
selection <- AICcmodavg::aictab(list.models)## Warning: no function found corresponding to methods exports from 'raster' for:
## 'wkt'selection #les deux premiers modèles semblent être ceux qui ont le plus de poids##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## buFishcond 5 392.62 0.00 0.74 0.74 -191.22
## buOnly 4 394.71 2.10 0.26 1.00 -193.30
## condOnly 4 408.99 16.38 0.00 1.00 -200.44
## null 3 412.61 20.00 0.00 1.00 -203.27
## NASConly 4 413.42 20.81 0.00 1.00 -202.65
## tempNASCFishcond 7 414.20 21.58 0.00 1.00 -199.93
## compONLY 4 414.48 21.86 0.00 1.00 -203.18
## tempONLY 4 414.61 22.00 0.00 1.00 -203.25est.bu <- AICcmodavg::modavg(list.models, "open_water_day")
est.bu##
## Multimodel inference on "open_water_day" based on AICc
##
## AICc table used to obtain model-averaged estimate:
##
## K AICc Delta_AICc AICcWt Estimate SE
## buFishcond 5 392.62 0.0 0.74 0.00 0
## buOnly 4 394.71 2.1 0.26 0.01 0
##
## Model-averaged estimate: 0
## Unconditional SE: 0
## 95% Unconditional confidence interval: 0, 0.01est.fish.cond <- AICcmodavg::modavg(list.models, "fish_cond")
est.fish.cond##
## Multimodel inference on "fish_cond" based on AICc
##
## AICc table used to obtain model-averaged estimate:
##
## K AICc Delta_AICc AICcWt Estimate SE
## tempNASCFishcond 7 414.20 21.58 0 0.18 0.08
## buFishcond 5 392.62 0.00 1 0.15 0.08
## condOnly 4 408.99 16.38 0 0.18 0.08
##
## Model-averaged estimate: 0.15
## Unconditional SE: 0.08
## 95% Unconditional confidence interval: 0.01, 0.3est.temp <- AICcmodavg::modavgShrink(list.models, "surf_temp_degC")## Warning in modavgShrink.AIClme(list.models, "surf_temp_degC"):
## Variables do not appear with same frequency across models, proceed with cautionest.temp##
## Multimodel inference on "surf_temp_degC" based on AICc
##
## AICc table used to obtain model-averaged estimate with shrinkage:
##
## K AICc Delta_AICc AICcWt Estimate SE
## null 3 412.61 20.00 0.00 0 0.00
## tempNASCFishcond 7 414.20 21.58 0.00 0 0.02
## buFishcond 5 392.62 0.00 0.74 0 0.00
## buOnly 4 394.71 2.10 0.26 0 0.00
## tempONLY 4 414.61 22.00 0.00 0 0.01
## condOnly 4 408.99 16.38 0.00 0 0.00
## NASConly 4 413.42 20.81 0.00 0 0.00
## compONLY 4 414.48 21.86 0.00 0 0.00
##
## Model-averaged estimate with shrinkage: 0
## Unconditional SE: 0
## 95% Unconditional confidence interval: 0, 0est.zoo <- AICcmodavg::modavgShrink(list.models, "NASC_zoo")## Warning in modavgShrink.AIClme(list.models, "NASC_zoo"):
## Variables do not appear with same frequency across models, proceed with cautionest.zoo##
## Multimodel inference on "NASC_zoo" based on AICc
##
## AICc table used to obtain model-averaged estimate with shrinkage:
##
## K AICc Delta_AICc AICcWt Estimate SE
## null 3 412.61 20.00 0.00 0 0
## tempNASCFishcond 7 414.20 21.58 0.00 0 0
## buFishcond 5 392.62 0.00 0.74 0 0
## buOnly 4 394.71 2.10 0.26 0 0
## tempONLY 4 414.61 22.00 0.00 0 0
## condOnly 4 408.99 16.38 0.00 0 0
## NASConly 4 413.42 20.81 0.00 0 0
## compONLY 4 414.48 21.86 0.00 0 0
##
## Model-averaged estimate with shrinkage: 0
## Unconditional SE: 0
## 95% Unconditional confidence interval: 0, 0new.data <- data.frame(
open_water_day = seq(from = min(dfish$open_water_day), to = max(dfish$open_water_day), length.out = 100),
fish_cond = mean(dfish$fish_cond),
region = "Amundsen Gulf Mouth",
surf_temp_degC = mean(dfish$surf_temp_degC),
NASC_zoo = mean(dfish$NASC_zoo),
indstandard = mean(dfish$indstandard)
)
preds <- AICcmodavg::modavgPred(cand.set = list.models, newdata = new.data)
preds2 <- predict(lme.glob2, newdata =new.data)
new.data$fit <- preds$mod.avg.pred
new.data$se <- preds$uncond.se
new.data$low95 <- preds$lower.CL
new.data$supp95 <- preds$upper.CL
#pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preds_bu.pdf")
plot(smooth(fit) ~ open_water_day, data = new.data,
ylab = "log10(Succès alimentaire des jeunes morues polaires)",
xlab = "Nombre de jours juliens depuis la débâcle",
main = "Valeurs prédites pondérées de succès alimentaire (+- IC 95%)",
type = 'l')
lines(x = new.data$open_water_day, y = smooth(new.data$supp95), lty = 'dotted')
lines(x = new.data$open_water_day, y = smooth(new.data$low95), lty = 'dotted')
#dev.off()Poster RSA QO
mean_logFS <- aggregate(logFS~ region + open_water_day, data = dfish, FUN = mean)
mean_logFS$region <- ordered(mean_logFS$region, c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "NEG"))
colours <- pals::stepped(n = 24)[c(1,4,9,11,12,13, 15,16,17)]
FS_OWD <- ggplot() +
geom_point(data = mean_logFS, mapping = aes(x = open_water_day, y = logFS, color = region), size = 10)+
geom_line(data = new.data, mapping = aes(y = fit, x = open_water_day)) +
geom_ribbon(data = new.data, aes(x = open_water_day, y = NULL, ymin = low95, ymax = supp95),
alpha = .15) +
theme_classic()+
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
legend.position = c(0.25,0.78),
axis.text.x = element_text(size=30),
axis.text.y = element_text(size = 30),
legend.text = element_text(size = 22),
axis.title = element_text(size = 30),
legend.title = element_text(size =22),
legend.background = element_rect(fill="transparent")) +
labs(x = "Nombre de jours depuis la débâcle", y = "log10(Succès alimentaire moyen)") +
scale_color_manual(values = colours, name = "Région", labels = c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "North-Eastern\nGreenland Sea")) +
scale_x_continuous(breaks = c(-20, 0, 20, 40, 60, 80, 100))
FS_OWD
ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/FS_OWD_fish.png", plot = FS_OWD, dpi = 500, width = 14, height = 10, bg = "transparent")library(pals)
colours_bar <- pals::stepped(n = 24)[c(20:17, 16:13, 12,5,8,4:1)]
gutcontent2$OWD_cat <- with(gutcontent2, ifelse(open_water_day <= 0, "before_bu",
ifelse(open_water_day >0 & open_water_day <20, "0-20",
ifelse(open_water_day >=20 & open_water_day<40, "20-40",
ifelse(open_water_day >=40 & open_water_day <60, "40-60",
ifelse(open_water_day >=60 & open_water_day < 80, "60-80", "80+"))))))
freqtable <- as.data.frame(table(gutcontent2$prey_id_comb,gutcontent2$region))
freqtable$Var2 <- ordered(freqtable$Var2, c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "NEG"))
freqtable$Var1 <- ordered(freqtable$Var1, c("appendicularia", "bivalve larvae", "egg", "limacina helicina", "calanoid sp c", "cyclopoidae sp c", "other copepodite sp", "other copepodite n", "other", "pseudocalanus sp c", "pseudocalanus sp n", "calanus glacialis c", "calanus hyperboreus c", "calanus sp c", "calanus sp n"))
preys_occ <- ggplot() +
geom_bar(data =freqtable, mapping = aes(fill=Var1, y=Freq, x=Var2), position="fill", stat="identity") +
ggtitle( "Fréquence d'occurence\ndans la diète")+
scale_fill_manual(name = "", values =colours_bar, labels = c("Appendicularia",
"Bivalvia",
"Oeuf",
"Gastropoda",
"Calanoidae spp.",
"Cyclopoidae spp.",
"Autre copépodite",
"Autre nauplii",
"Autre type de proies",
expression(paste(italic("Pseudocalanus "), "spp.", sep = "")),
expression(paste(italic("Pseudocalanus "), "spp. Nauplii", sep = "")),
expression(italic("Calanus glacialis")),
expression(italic("Calanus hyperboreus")),
expression(paste(italic("Calanus "), "spp.", sep = "")),
expression(paste(italic("Calanus "), "spp. Nauplii", sep = ""))
))+
theme_classic() +
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size = 20),
legend.text = element_text(size = 22),
legend.text.align = 0,
axis.title = element_text(size = 25),
legend.title = element_text(size =22),
legend.background = element_rect(fill=alpha('white', 0.4)),
plot.title = element_text(size = 30, hjust = 0.5)) +
labs(y = "Proportion", x = "") +
scale_x_discrete(labels = c("MS", "AGM", "CM", "LV", "PS", "LS", "NW", "WBB", "NEG"))
preys_occ
#ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preys_occ.png", plot = preys_occ, dpi = 500, width = 14, height = 10, bg = "transparent")
carbon_input <- aggregate(carbon_mg ~ prey_id_comb + region, data = gutcontent2, FUN = sum)
carbon_input$region <- ordered(carbon_input$region, c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "NEG"))
carbon_input$prey_id_comb <- ordered(carbon_input$prey_id_comb, c("appendicularia", "bivalve larvae", "egg", "limacina helicina", "calanoid sp c", "cyclopoidae sp c", "other copepodite sp", "other copepodite n", "other", "pseudocalanus sp c", "pseudocalanus sp n", "calanus glacialis c", "calanus hyperboreus c", "calanus sp c", "calanus sp n"))
preys_carbon <- ggplot() +
geom_bar(data =carbon_input, mapping = aes(fill=prey_id_comb, y=carbon_mg, x=region), position="fill", stat="identity") +
ggtitle("Proportion en carbone\ndans la diète")+
scale_fill_manual(name = "", values =colours_bar, labels = c("Appendicularia",
"Bivalvia",
"Oeuf",
"Gastropoda",
"Calanoidae spp.",
"Cyclopoidae spp.",
"Autre copépodite",
"Autre nauplii",
"Autre type de proies",
expression(paste(italic("Pseudocalanus "), "spp.", sep = "")),
expression(paste(italic("Pseudocalanus "), "spp. Nauplii", sep = "")),
expression(italic("Calanus glacialis")),
expression(italic("Calanus hyperboreus")),
expression(paste(italic("Calanus "), "spp.", sep = "")),
expression(paste(italic("Calanus "), "spp. Nauplii", sep = ""))
)) +
theme_classic() +
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size = 20),
legend.text = element_text(size = 22),
legend.text.align = 0,
axis.title = element_text(size = 25),
legend.title = element_text(size =22),
legend.background = element_rect(fill=alpha('white', 0.4)),
plot.title = element_text(size = 30, hjust = 0.5)) +
labs(y = "", x = "Région") +
scale_x_discrete(labels = c("MS", "AGM", "CM", "LV", "PS", "LS", "NW", "WBB", "NEG"))
preys_carbon
#ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preys_carbon.png", plot = preys_carbon, dpi = 500, width = 14, height = 10, bg = "transparent")
freqtable3 <- as.data.frame(table(gutcontent2$prey_id_comb,gutcontent2$OWD_cat))
freqtable3$Var1 <- ordered(freqtable3$Var1, c("appendicularia", "bivalve larvae", "egg", "limacina helicina", "calanoid sp c", "cyclopoidae sp c", "other copepodite sp", "other copepodite n", "other", "pseudocalanus sp c", "pseudocalanus sp n", "calanus glacialis c", "calanus hyperboreus c", "calanus sp c", "calanus sp n"))
freqtable3$Var2 <- ordered(freqtable3$Var2,c("before_bu", "0-20","20-40", "40-60", "60-80", "80+"))
preys_occ_owd <- ggplot() +
geom_bar(data =freqtable3, mapping = aes(fill=Var1, y=Freq, x=Var2), position="fill", stat="identity") +
ggtitle( "Fréquence d'occurence\ndes proies dans la diète")+
scale_fill_manual(name = "", values =colours_bar, labels = c("Appendicularia",
"Bivalvia",
"Oeuf",
"Gastropoda",
"Calanoidae spp.",
"Cyclopoidae spp.",
"Autre copépodite",
"Autre nauplii",
"Autre type de proies",
expression(paste(italic("Pseudocalanus "), "spp.", sep = "")),
expression(paste(italic("Pseudocalanus "), "spp. Nauplii", sep = "")),
expression(italic("Calanus glacialis")),
expression(italic("Calanus hyperboreus")),
expression(paste(italic("Calanus "), "spp.", sep = "")),
expression(paste(italic("Calanus "), "spp. Nauplii", sep = ""))
))+
theme_classic() +
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
axis.text.x = element_text(size=30),
axis.text.y = element_text(size = 30),
legend.text = element_text(size = 22),
legend.text.align = 0,
axis.title = element_text(size = 35),
legend.title = element_text(size =22),
legend.background = element_rect(fill=alpha('white', 0.4)),
plot.title = element_text(size = 35, hjust = 0.5)) +
labs(y = "Proportion", x = "Nombre de jours depuis la débâcle")+
scale_x_discrete(labels = c("Avant\ndébâcle", "0-20 j.", "20-40 j.", "40-60 j.", "60-80 j.", "80+ j."))
preys_occ_owd
ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preys_OWD.png", plot = preys_occ_owd, dpi = 500, width = 17, height = 14, bg = "transparent")
carbon_input2 <- aggregate(carbon_mg ~ prey_id_comb + OWD_cat + unique_fish_id, data = gutcontent2, FUN = sum)
carbon_input3 <- aggregate(carbon_mg ~ prey_id_comb + OWD_cat, data = carbon_input2, FUN=mean)
carbon_input3$OWD_cat <- ordered(carbon_input3$OWD_cat ,c("before_bu", "0-20","20-40", "40-60", "60-80", "80+"))
carbon_input3$prey_id_comb <- ordered(carbon_input3$prey_id_comb, c("appendicularia", "bivalve larvae", "egg", "limacina helicina", "calanoid sp c", "cyclopoidae sp c", "other copepodite sp", "other copepodite n", "other", "pseudocalanus sp c", "pseudocalanus sp n", "calanus glacialis c", "calanus hyperboreus c", "calanus sp c", "calanus sp n"))
preys_carbon_OWD <- ggplot() +
geom_bar(data =carbon_input3, mapping = aes(y=carbon_mg, x=OWD_cat, fill = prey_id_comb), stat="identity") +
ggtitle("Moyenne de carbone ingéré par poisson\net contribution des proies")+
scale_fill_manual(name = "", values =colours_bar, labels = c("Appendicularia",
"Bivalvia",
"Oeuf",
"Gastropoda",
"Calanoidae spp.",
"Cyclopoidae spp.",
"Autre copépodite",
"Autre nauplii",
"Autre type de proies",
expression(paste(italic("Pseudocalanus "), "spp.", sep = "")),
expression(paste(italic("Pseudocalanus "), "spp. Nauplii", sep = "")),
expression(italic("Calanus glacialis")),
expression(italic("Calanus hyperboreus")),
expression(paste(italic("Calanus "), "spp.", sep = "")),
expression(paste(italic("Calanus "), "spp. Nauplii", sep = ""))
))+
theme_classic() +
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
axis.text.x = element_text(size=28),
axis.text.y = element_text(size = 28),
legend.text = element_text(size = 25),
legend.text.align = 0,
axis.title = element_text(size = 35),
legend.title = element_text(size =22),
legend.background = element_rect(fill=alpha('white', 0.4)),
plot.title = element_text(size = 35, hjust = 0.5))+
scale_x_discrete(labels = c("Avant\ndébâcle", "0-20 j.", "20-40 j.", "40-60 j.", "60-80 j.", "80+"))+
labs(y = "Moyenne de carbone ingéré/poisson (mg)", x = "Nombre de jours depuis la débâcle")
preys_carbon_OWD
ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preys_carbon_OWD.png", plot = preys_carbon_OWD, dpi = 500, width = 17, height = 14, bg = "transparent")
#le faire pour des catégories de succès alimentaire
carbon_input4 <- aggregate(carbon_mg ~ prey_id_comb + open_water_day + unique_fish_id, data = gutcontent2, FUN = sum)
carbon_input5 <- aggregate(carbon_mg ~ prey_id_comb + open_water_day, data = carbon_input4, FUN = mean)
carbon_input6 <- aggregate(carbon_mg ~ open_water_day, data = carbon_input5, FUN = sum)
carbon_input_caspC <- aggregate(carbon_mg ~ open_water_day,
data = carbon_input4[carbon_input5$prey_id_comb %in%
c("calanus glacialis c",
"calanus hyperboreus c",
"calanus sp c"),], FUN = sum)
ggplot()+
geom_point(data = carbon_input6, mapping = aes(x=open_water_day, y = carbon_mg), color = "black")+
geom_point(data = carbon_input_caspC, mapping = aes(x=open_water_day, y = carbon_mg), color = "#990F26")+
geom_smooth(method = "loess", span = 0.75,data = carbon_input6, mapping = aes(x=open_water_day, y = carbon_mg), color = "black", fill = "grey") +
geom_smooth(method = "loess", span = 0.75, data = carbon_input_caspC, mapping = aes(x=open_water_day, y = carbon_mg), color = "#990F26", fill = "#990F26")+
theme_classic() +
theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA),
axis.text.x = element_text(size=28),
axis.text.y = element_text(size = 28),
legend.text = element_text(size = 22),
legend.text.align = 0,
axis.title = element_text(size = 30),
legend.title = element_text(size =22),
legend.background = element_rect(fill=alpha('white', 0.4)),
plot.title = element_text(size = 30, hjust = 0.5))+
scale_x_continuous(breaks = c(-20, 0, 20, 40, 60, 80, 100))+
scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1,1.25))